///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/fqmc/accumulator.h
///
/// \brief   header file for the Accumulator base class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef FORMIC_FQMC_ACCUMULATOR_HEADER
#define FORMIC_FQMC_ACCUMULATOR_HEADER

#include<string>

#include<boost/shared_ptr.hpp>

#include<formic/mpi/interface.h>
#include<formic/wfn/coeff_obj.h>
#include<formic/numeric/numeric.h>

namespace formic {

  class InputBase;
  class Configuration;

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  the base class for accumulator objects that collect data during a monte carlo sample
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class Accumulator {

    protected:

      /// \brief  the coefficient object of the guiding wave function
      boost::shared_ptr<formic::CoeffObj<S> > _base_gf;

      /// \brief  the coefficient object of the value wave function
      boost::shared_ptr<formic::CoeffObj<S> > _base_vf;

      /// \brief  total weight accumulated at all configurations
      double _total_weight;
      double _naccept;

      /// \brief  average square value function to guiding function ratio
      double _vgsa;

      /// \brief  local value of the square of the value function to guiding function ratio
      double _vgs;

      /// \brief  binary operator function object used in computing the variance
      struct VarBinOp {

        double m;

        VarBinOp(double mean) : m(mean) {}
      
        double operator()(double x, double y) { return x + y*y - m*m; }
        
      };

      /// \brief  binary operator function object used in computing the autocorrelation
      struct AutoCorrBinOp {

        double m;
        int t;

        AutoCorrBinOp(double mean, int offset_time) : m(mean), t(offset_time) {}

        double operator()(const double x, const double & y) {
          return x + ( y - m ) * ( *(&y-t) - m );
        }

      };

    public:

      /// \brief  default constructor
      Accumulator() : _total_weight(0.0), _naccept(0.0), _vgsa(0.0), _vgs(0.0) {}

      /// \brief  constructor with guiding and value functions
      Accumulator(const boost::shared_ptr<formic::CoeffObj<S> > & gf,
                  const boost::shared_ptr<formic::CoeffObj<S> > & vf)
        : _total_weight(0.0),
          _naccept(0.0),
          _vgsa(0.0),
          _vgs(0.0),
          _base_gf(gf),
          _base_vf(vf)
      {}

      /// \brief  make destructor virtual so that deleting a pointer-to-base will call childrens' destructors
      virtual ~Accumulator() {}

      virtual std::string name() const = 0;

      /// \brief  return the acceptance ratio
      double acceptance_ratio() const { assert( _total_weight != 0.0 ); return _naccept / _total_weight; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Accumulate data for the given configuration.
      ///
      /// \param[in]      userinp  the user's input options
      /// \param[in]      w        the weight to apply to this configuration's contribution
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void accumulate(const formic::InputBase & userinp, const double w) {

        // ensure weight is positive
        if ( w < 0.0 )
          throw formic::Exception("unexpected negative weight in %s") % this->name();

        // add the weight of the configuration to the total
        _total_weight += w;

        // count the acceptance
        _naccept += 1.0;

        // get the square value to guiding ratio and add it to the total
        _vgs = ( _base_vf == _base_gf  ?  1.0  :  formic::square_norm( ( _base_vf->get_coeff() / _base_gf->get_coeff() ).get() ) );
        _vgsa += w * _vgs;

        // perform any special accumulations as specified by the child class
        this->child_accumulate(userinp, w);

      }

      virtual void child_accumulate(const formic::InputBase & userinp, const double w) = 0;

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief finalizes the accumulation (this may include for example communicating the results)
      ///
      /// \param[in]      userinp  the user's input options
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void finalize(const formic::InputBase & userinp) {

        // set up space for mpi reduction
        const int nred = 3;
        double x[nred];
        double y[nred];

        // prepare for reduce
        x[0] = _total_weight;
        x[1] = _vgsa;
        x[2] = _naccept;

        // all reduce
        formic::mpi::allreduce(x, y, nred, MPI::SUM);

        // record results
        _total_weight = y[0];
        _vgsa = y[1] / _total_weight;
        _naccept = y[2];

        // perform any finalizations required by the child
        this->child_finalize(userinp);

      }

      virtual void child_finalize(const formic::InputBase & userinp) = 0;

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  resets the accumulator
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void reset(const formic::Configuration & config) {

        // reset the total weight, the average square value to guiding ratio, and number of accepts
        _total_weight = 0.0;
        _vgsa = 0.0;
        _naccept = 0.0;

        // perform any resetting required by the child
        this->child_reset(config);

      }

      virtual void child_reset(const formic::Configuration & config) = 0;

  };

}

#endif
